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Abstract 

Transient  state  response  analysis  of  phosphoric  acid  fuel  cell  (PAFC)  cathode  is  important  to  understand  various  competitive  processes 
like  diffusion,  reaction  and  product  back  diffusion  occurring  at  various  layers  of  the  composite  cathode.  A  two-dimensional  unsteady  state 
model  for  simulating  PAFC  cathode  is  developed  as  an  extension  of  the  previously  developed  steady-state  model  [S.  Roy  Choudhury,  M.B. 
Deshmukh,  R.  Rengaswamy,  A  two-dimensional  steady-state  model  for  phosphoric  acid  fuel  cells  (PAFC),  J.  Power  sources  112  (2002) 
137-152].  The  transient  model  is  solved  to  study  the  impact  of  various  parameters  such  as  Tafel  slope,  diffusivity  etc  on  the  step  response  of 
the  fuel  cell.  The  effect  of  partial  pressure  variation  in  bulk  gas  for  large  sized  PAFC  cathode  is  also  analysed.  Trend  analysis  based  on  the 
model  output  is  also  experimentally  verified  using  a  small  unit  cell  setup.  The  effect  of  various  parameters  on  the  settling  time  of  the  cathode, 
as  revealed  in  this  study,  suggests  possible  development  of  a  diagnostic  tool  employing  such  transient  model. 

©  2004  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Owing  to  diverse  applications  of  phosphoric  acid  fuel  cells 
(PAFCs),  dynamic  operation  of  a  fuel  cell  stack  is  quite  com¬ 
mon.  Such  applications  includes  vehicular  propulsion,  hybrid 
backup  systems,  several  military  applications  etc,  where  the 
load  keeps  fluctuating  unlike  steady  on-site  power  generation 
systems. 

A  PAFC  is  composed  of  two  porous  gas  diffusion  elec¬ 
trodes,  the  anode  and  cathode.  The  gas  after  passing  through 
porous  substrate  (diffusion  layer)  reaches  the  catalyst  layer 
also  termed  as  reaction  layer.  The  catalyst  layer  is  a  porous 
layer,  partly  filled  with  electrolyte.  The  gas  diffuses  through 
the  dry  pores,  dissolves  in  the  electrolyte  and  then  diffuses  to 
the  catalytic  site.  Oxygen  reduction  takes  place  on  the  site  and 
the  product  water  back  diffuses  through  the  diffusion  layer 
to  escape  in  the  oxygen  stream. 


*  Corresponding  author.  Tel.:  +1  315  268  4423;  fax:  +1  315  268  6654. 
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Unsteady  state  simulation  of  fuel  cell  is  necessary  to  un¬ 
derstand  the  system  transient  response  under  variable  loads 
with  various  time  constants.  Apart  from  that,  unsteady  state 
simulation  model  is  also  important,  to  analyze  the  effect 
of  disturbances  of  input  variables  even  for  a  steady- state 
operation.  A  good  representative  transient  state  model  can 
be  used  for  diagnostics  of  a  running  stack,  which  is  very 
important  for  high  reliability  applications  such  as  hospital 
and  military  applications.  The  unsteady  state  model  further 
finds  its  use  in  analyzing  several  control  strategies  to  operate 
fuel  cells,  in  the  best  possible  way  under  diverse  operating 
conditions. 


2.  Development  of  unsteady  state  model 

The  steady-state  model  for  a  PAFC  (Fig.  1)  presented  in 
Choudhury  et  al.  [1]  is  further  extended  to  develop  an  un¬ 
steady  state  model.  Bulk  gas,  i.e.,  oxygen  flows  in  the  In¬ 
direction  through  the  ribbed  support  over  the  left  side  of 
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Nomenclature 
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ac 
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cross  section  area  of  an  agglomerate  in  YZ 
plane  [cm2] 

cross  section  area  of  each  groove  [cm2] 
effective  area  of  catalyst  site  per  unit  volume 
of  agglomerate  [cm2  cm-3] 
groove  cross  section  [sqcm] 
normal  Tafel  slope  [V  per  decade] 
concentration  of  oxygen  on  surface  of  agglom¬ 
erate  [gmol  cm3  ] 

concentration  of  oxygen  inside  agglomerate 
[gmol  cm-3] 

phosphoric  acid  concentration  in  percent  ortho 
phosphoric  acid  equivalent 
diffusivity  inside  agglomerate  for  oxygen 
[cm2  s-1] 

diffusivity  inside  liquid  encapsulation  on  ag¬ 
glomerate  for  oxygen  [cm2  s-1] 
diffusivity  in  diffusion  layer  for  oxygen- water 
[cm2  s-1] 

diffusivity  in  reaction  layer  for  oxygen-water 
[cm2  s-1] 

hydraulic  dia  of  each  groove  [cm] 
open  circuit  potential  [V] 
standard  open  circuit  potential  [V] 
oxygen  electrode  (solid  carbon  phase)  poten¬ 
tial  [V] 

local  electrolyte  (liquid)  potential  [V] 
friction  factor 

Faraday  constant  [96,500  coulomb  g-1  equiv¬ 
alent  weight] 

Henry’s  constant  for  oxygen  solubility  in  phos¬ 
phoric  acid  [gmol  cm-3  atm-1] 
local  current  density  [A  cm-2] 
exchange  current  density  for  oxygen  reduction 
on  Pt  [A  cm-2] 
current  in  Y- direction  [A] 
flow  of  cation  in  Y-direction  [A] 
flow  of  cation  in  F-direction  [A] 
reaction  constant  of  oxygen  reduction  reaction 
groove  length  [cm] 

product  of  molar  flux  of  ith  species  at  gas  elec¬ 
trode  interphase  in  Y-direction  and  width  of 
the  cell  [mol  cm-1  s-1] 
inlet  molar  flow  rate  of  ith  species  [gmol  s-1] 
number  of  electrons  taking  part  in  the  reaction 
grooves  per  manifold 

number  of  agglomerate  exposed  per  unit  sqcm 
area  in  YZ  plane 

oxygen  reduction  rate  [gmol  s-1  cm-3] 
total  pressure  (atm) 

partial  pressure  of  the  ith  species  (atm) 


Pit) 

partial  pressure  of  the  ith  species  at  inlet  (atm) 

R 

universal  gas  constant  [atm  cm3  gmol-1  K-1] 

^agg 

average  radius  of  an  agglomerate  [cm] 

r 

radial  dimension  inside  an  agglomerate  origi¬ 
nating  at  the  center  [cm] 

T 

system  temperature  in  Kelvin  scale 

t 

time  [s] 

U 

oxygen  reaction  potential  w.r.t.  the  reference 
[V] 

W 

width  of  the  cell  in  Z-direction  (cm) 

fF(max) 

maximum  electrical  work  available  [W] 

/SE° 

^^rev 

standard  half  cell  potential  w.r.t.  hydrogen 
electrode  [V] 

a 

activation  coefficient 

€ 

capacitance  per  unit  volume  of  reaction  layer 
[farad  cm-3] 

0 

local  overpotential  [V] 

00 

overpotential  at  the  bulk-electrolyte  and  elec¬ 
trode  interphase  [V] 

Ke 

electrolyte  conductivity  [mho  cm-1] 

< 

electrolyte  conductivity  inside  agglomerate 
core  [mho  cm-1] 

A-d 

porosity  of  diffusion  layer 

xr 

porosity  of  reaction  layer 

Pa. 

agglomerate  density  inside  reaction  layer 
[number  cm-3] 

V 

flow  per  manifold  at  system  temperature  and 
pressure  [cc  s-1] 

Subscripts  for  ith  species 

1 

oxygen 

2 

water 

the  electrode.  As  the  gas  flows,  part  of  the  oxygen  diffuses 
through  the  porous  substrate  of  the  electrode.  The  diffusion 
of  oxygen  is  in  Y-direction. 

Diffusion  of  oxygen  and  back  diffusion  of  water  vapor  in 
the  porous  substrate  for  unsteady  state  can  thus  be  defined  by 
equation  1  which  is  as  follows 


Similarly,  inside  the  reaction  layer,  mass  transfer  with  re¬ 
action  can  be  described  by  Eq.  (2) 


(2) 
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where 

r  =  47T/9aZ)agg7?aggCos[0(7?agg)  coth  4>(Ra gg)  —  1] 

0(^agg)  = 

and 

cos  =  Hpi 


(3) 

(4) 

(5) 


It  is  assumed  that  the  agglomerates  being  very  small  in  size 
will  react  very  fast  to  the  local  surface  gas  partial  pressure 
and  hence  the  unsteady  state  contribution  of  the  agglomerates 
are  not  considered. 

Distribution  of  overpotential  inside  the  catalyst  layer  can 
be  modeled  as  done  for  the  steady-state  model,  by  considering 
ohmic  law  for  the  cation  flow  [2],  in  this  case  hydrogen  ions, 
in  the  electrolyte  filled  zone.  Thus,  change  of  overpotential 
in  X-  and  T-directions  can  be  related  to  the  current  flux  as 
follows 

d2i]  d2t]  ( nFr )  f  dr]\ 

— o  +  — T  —  -  T  (  —  1  (6) 

dx 2  dy 2  KqA*N1  \  3 1 ) 

The  capacitance  per  unit  volume  of  catalyst  layer  aac 
arises  primarily  from  the  double  layer  near  the  active  Pt 
sites. 


3.  Boundary  conditions 

Providing  boundary  conditions  for  solving  the  unsteady 
state  simulation  depends  upon  the  mode  of  operation,  starting 
conditions  etc. 

Following  are  two  most  common  mode  of  operation: 


Type-1 :  The  system  is  running  under  steady  state,  and  from 
that  one  or  more  operating  parameters  like  gas  flowrate, 
electrical  load,  humidity  level  is  changed  with  respect 
to  time. 

Type-2 :  Inlet  and  outlet  of  the  system  is  closed  and  heated 
up  to  the  required  temperature  and  then  suddenly  the 
valves  were  open  to  start  the  cell  with  operating  pa¬ 
rameters  varied  with  respect  to  time. 

Type- 1  operation  is  of  primary  importance  as  it  is  the  con¬ 
dition  that  is  valid  once  the  cell  has  started.  Information  re¬ 
garding  the  initial  values  of  the  dependent  variables  over  the 
region  at  time  t  —  0  is  provided  from  the  steady- state  simu¬ 
lation  of  the  steady  state  from  which  the  changeover  of  op¬ 
erating  parameters  is  considered. 

Boundary  conditions  numbered  1-10  as  mentioned  in 
Choudhury  et  al.  [1]  for  time  >  0  is  used  as  well.  Depend¬ 
ing  upon  the  nature  of  variation  of  operating  conditions  with 
respect  to  time,  the  particular  conditions  are  replaced  with 
the  time  dependent  profile  of  the  same.  For  e.g.,  if  the  load  is 
changed  suddenly,  then  condition  no.  7  [1]  may  be  modified 
as  follows: 

rj  =  r)o  at  t  <  0  and  ij  =  r]i  at  t  >  0  at  v  =  B  Wy 

In  case  of  change  of  inlet  flowrate  or  concentration,  the  par¬ 
ticular  boundary  conditions  are  modified  accordingly. 

3.1.  Experimental  setup  for  dynamic  data  acquisition 

The  experimental  setup  for  validation  of  the  unsteady  state 
model  is  similar  to  that  of  the  unit  cell  study  as  mentioned  in 
Choudhury  et  al.  [1].  A  computerised  data  acquisition  sys¬ 
tem  is  added  to  record  the  transient  output  of  the  cell.  The 
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Fig.  2.  Unit  cell  experimentation  setup  for  transient  response  measurement. 


compensator  circuit  is  replaced  with  a  high  current  electro¬ 
chemical  workstation  (AutoLAB) which  works  as  a  variable 
load  and  any  potential  profile  can  be  generated  with  this  sys¬ 
tem.  A  secondary  PC  monitors  the  cell  performance.  Fig.  2 
shows  the  complete  experimental  setup.  A  unit  cell  setup  was 
used  as  the  ionic  transient  potential  is  required  to  be  mea¬ 
sured  to  check  the  capacitance  effect.  The  data  acquisition 
card  used  is  analog  subsystem,  FIFO,  non-multiplexing  type, 
from  Advantech,  USA.  It  has  a  maximum  sampling  speed  of 
100  sample  per  second  per  channel. 


4.  Analysis  of  the  simulated  results 

The  unsteady  state  model  is  solved  using  the  software, 
Pdease2D  that  is  used  for  solving  the  steady- state  model 
[1].  The  time  dependent  part  is  solved  by  marching  forward 
method.  After  a  small  time  increment,  the  partial  pressure 
along  the  bulk  gas  flow  interphase  is  updated  as  mentioned 
in  the  solution  to  the  steady  state  model  [1].  In  this  work,  the 
effect  of  variations  in  the  model  parameters  on  step  response 
of  the  electrode  is  tested.  Values  of  the  parameters  assumed 
are  provided  in  the  Table  1,  unless  mentioned  otherwise  in 
the  figures. 

The  simulated  response  is  normalized  by  the  final  steady- 
state  value.  This  is  done  in  order  to  understand  the  response 
time  better  by  making  the  ordinate  scale  same.  We  first  report 
the  results  for  the  one-dimensional  model.  Fig.  3  depicts  the 
simulated  one-dimensional  current  response  over  a  step  jump 
from  0. 1  V  to  0.2  V  overpotential  at  varying  tafel  slopes.  The 
figure  shows  that  as  the  tafel  slope  decreases,  the  peak  in¬ 
creases,  although  the  curve  becomes  more  steeper,  and  the 
time  required  to  attain  steady- state  value  also  decreases.  In¬ 
crease  of  the  peak  current  with  decreased  tafel  slope  is  ex¬ 
pected,  as  sudden  increase  in  overpotential  along  with  the 
local  high  reactant  concentration  generates  higher  current. 
This  effect  also  consumes  the  excess  reactant  faster  and  sets 
up  the  diffusion  pattern  early  so  that  the  settling  time  is  de¬ 


creased.  Thus,  the  model  result  matches  with  the  expected 
behavior. 

Fig.  4  depicts  the  simulated  step  response  for  two  val¬ 
ues  of  reaction  layer  diffusivity.  As  the  diffusion  coefficient 
of  the  reaction  layer  increases,  the  normalized  peak  current 
increases  along  with  the  settling  time.  Due  to  enhanced  dif¬ 
fusion  properties,  the  steady-state  current  is  higher.  Thus,  the 
normalized  current  which  is  a  ratio  between  the  transient  cur¬ 
rent  and  the  final  steady-state  current  shows  a  smaller  peak  at 
the  potential  changeover  time  for  the  lower  diffusivity  case. 
Owing  to  better  diffusion  characteristics,  the  settling  time  is 

Table  1 


Values  of  parameters 


Parameter  name 

Value 

Remark 

A* 

7r/?2gg  cm2  of  projection  in  YZ 
plane 

Literature 

<2a 

3.5E04cm2  cm-3  of  agglom¬ 
erate 

Experimental 

b 

60  mV  per  decade 

Tuning 

^agg 

1.0E— 04  cm2  s-1 

Literature 

Dr 

0.2E-03  cm2  s_1 

Literature 

Diffusion  layer 

400  [xm  thick 

Experimental 

E 

0.92  V 

Experimental 

H 

5.5E— 08  gmolcm-3  atm-1 

Literature 

io 

1.0E— 06  A  cm-2 

Literature 

n 

4 

Stoichiometry 

VI 

0.5//7«2gg 

Literature  [2] 

^agg 

1.0E— 04  cm 

SEM  picture 

Reaction  layer 

400  [xm  thick 

Experimental 

T 

150  +  273  K 

Experimental 

W 

10  cm 

Experimental 

Kt 

Function  of  phosphoric  acid 
concentration 

Model 

Pa 

(3/4)[0.5/(jri?agg)] 

Literature 

A-d 

0.5 

Experimental 

Ar 

0.5 

Experimental 

ac 

0.2  sqcm 

Experimental 

ns 

13  grooves  per  manifold 
(Type-A) 

Experimental 

^agg 

Dagg/Ar 

Literature 

< 

0.25 

Literature 
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Fig.  3.  Simulated  step  response  at  different  tafel  slope  ( b  =  60  and  90 
mV/decade)  in  one-dimensional  mode. 


Fig.  6.  Simulated  step  response  in  two-dimensional  mode. 
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Fig.  4.  Simulated  step  response  at  different  diffusivity  value  in  one¬ 
dimensional  mode. 

increased  as  better  “supply”  of  reactant  allows  slow  reduction 
of  the  current. 

After  checking  the  one-dimensional  model,  the  two- 
dimensional  model  is  simulated.  Fig.  5  depicts  the  steady- 
state  oxygen  profile  prior  to  the  step  change.  Here,  a  simple 
"Type- A"  groove,  as  mentioned  in  [1],  which  is  composed 


Fig.  5.  Simulated  oxygen  profile  prior  to  a  step  jump  in  two-dimensional 
mode. 


Fig.  7.  Simulated  one-  and  two-dimensional  step  response. 

of  all  parallel  channels  is  used.  This  groove  has  already  been 
described  in  the  steady-state  model  validation  section  of  [1]. 
The  flow  rate  used  is  1  lm-1  of  oxygen  flow  per  cell.  Fig.  6 
depicts  the  step  response  of  the  model.  This  response  if  plot¬ 
ted  with  the  same  parameters  but  using  the  one-dimensional 
model  shows  that  the  behavior  is  virtually  the  same  as  shown 
in  the  Fig.  7.  This  feature  could  be  quite  helpful  for  diagnostic 
purposes.  For  example,  if  the  DC  polarization  performance 
of  a  large  format  cell  decreases,  and  there  is  also  a  change 
in  settling  time,  it  implies  that  the  possible  reason  could  be 
some  problem  in  the  electrode  rather  than  partial  blocking  or 
low  gas  flow  rate. 

5.  Experimental  results 

5.7.  Varying  reaction  layer  diffusivity 

Reaction  layer  diffusion  is  highly  important  as  the  PAFC 
cathode  performance  depends  largely  on  the  supply  of  reac¬ 
tants  to  the  active  sites  and  back  diffusion  of  water.  As  the 
cathode  ages,  owing  to  mild  carbon  corrosion  at  the  teflon- 
carbon  interphase,  more  active  area  is  wetted  by  the  phos¬ 
phoric  acid.  This  increases  the  total  active  area  inside  the  re¬ 
action  layer.  However,  as  the  wetting  increases,  the  dry  pores 
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Fig.  9.  Experimental  step  response  with  simulated  data. 


Fig.  8.  Experimental  step  response  of  unit  cell  cathode  of  0.1  V  to  0.2  V 
overpotential  change. 

inside  the  reaction  layer  decreases.  This  results  in  lowering 
of  effective  gas  diffusivity.  Initially,  as  the  cathode  ages,  the 
performance  increases  although  the  diffusivity  is  lowered. 
However,  as  the  wetting  passes  a  certain  limit,  due  to  severe 
drop  in  diffusivity,  the  overall  cathode  performance  start  de¬ 
creasing. 

The  effect  of  cathode  reaction  layer  diffusivity,  on  settling 
time,  is  studied  using  the  unit  cell  setup  as  described  before. 
The  cathode  is  prepared  as  mentioned  in  [1].  The  cell  after 
soaking  in  88%  phosphoric  acid  for  24  h  at  a  temperature  of 
150  °C,  is  assembled  and  operated  for  4  h  at  a  current  density 
of  200  mA  cm-2.  All  the  readings  are  taken  after  this  initial 
conditioning  of  the  cathode  to  stabilize  the  system.  In  order 
to  age  the  cathode,  so  that  the  diffusivity  of  the  electrode  re¬ 
duces,  a  rapid  aging  method  is  implemented.  The  cell  is  kept 
at  open  circuit  mode  for  several  hours  after  the  initial  con¬ 
ditioning.  The  step  response  analysis  is  done  after  different 
hours  of  aging  to  understand  the  trend. 

Fig.  8  depicts  the  result  of  the  step  response.  It  can  be 
seen  that  as  the  cathode  ages,  the  settling  time  decreases. 
This  results  matches  quite  well  with  the  simulation  result  as 
depicted  in  Fig.  4. 

5.2.  Matching  of  electrode  response  with  simulated  data 

The  normalized  experimental  result  of  the  unit  cell  as  men¬ 
tioned  in  the  Section  3.1  is  plotted  along  with  various  other 
simulated  data  in  the  Fig.  9.  It  can  be  seen  that  the  result  is 
matching  when  the  tafel  slope  is  60  mV  per  decade,  which 
is  that  of  the  polished  Pt  surface,  with  the  other  parameters 
remaining  same  as  in  the  Table  1.  However,  the  previously 
tuned  data  for  steady-state  analysis  [1]  shows  a  tafel  slope  of 
around  120  mV  per  decade.  This  discrepancy  could  be  ana¬ 


lyzed  as  follows.  The  tafel  slope  of  a  rough  catalyst  surface 
as  expected  in  the  PAFC  cathode  catalyst  is  higher  than  the 
polished  Pt  surface  due  to  near  surface  convection  which  is 
not  modeled.  Only  bulk  diffusion  inside  the  agglomerate  is 
considered  while  developing  the  agglomerate  model.  For  the 
short  transient  period,  the  response  thus  matches  with  the  real 
tafel  slope  value. 


6.  Conclusion 

The  models  developed  are  studied  for  step  response.  It  is 
observed  that  the  settling  time  depends  on  several  electrode 
related  parameters  such  as  diffusivity,  tafel  slope  etc.  How¬ 
ever,  the  settling  time  remains  the  same  for  concentration 
gradient  along  the  bulk  flow,  which  is  external  to  the  elec¬ 
trode.  This  aspect  may  be  valuable  for  diagnostic  purposes. 
For  example,  same  settling  time  of  a  normal  cell  and  a  faulty 
cell  might  indicate  an  external  reactant  blockage  in  the  gas 
flow  grooves,  and  the  basic  cell  is  not  affected.  As  a  cell  ages, 
the  wettability  of  the  electrode  increases,  thus  increasing  the 
active  area  of  the  electrode,  but  decreasing  the  reaction  layer 
diffusivity.  The  extent  of  electrode  wetting  in  a  newly  as¬ 
sembled  cell  can  be  ascertained  though  the  step  response,  as 
suggested  by  the  simulation  in  Fig.  4. 
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